COVID-19 Progression

library(tidyverse)
library(magrittr)
library(lubridate)
library(scales)
library(matrixStats)
library(ggrepel)
library(broom)
library(glue)
library(jsonlite)
library(rvest)
library(RCurl)
library(pander)
library(DT)
library(plotly)
panderOptions("big.mark", ",")
panderOptions("table.split.table", Inf)
panderOptions("table.style", "rmarkdown")
panderOptions("missing", "")
theme_set(theme_bw())

# Handle updates between 12am & 12pm
dt <- Sys.Date()
if (as.numeric(format(Sys.time(), "%H")) <= 12){
  dt <- Sys.Date() - 1
}

Disclaimer: This very simple report was prepared by a bioinformatician with no experience in epidemiology or virology, and as such should be treated simply as an alternate viewpoint on the data, which I was simply unable to find elsewhere. Many other people exist with much greater expertise on this subject. However, I do hope this provides a useful perspective which is able to add constructively to the wider discussion. In addition, it should be noted that this is very much focussed on Australian data.

Data Sources

confirmed <- url("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv") %>%
    read_csv() %>%
    pivot_longer(
        cols = ends_with("20"),
        names_to = "date",
        values_to = "confirmed"
    )  %>%
    mutate(
        date = str_replace_all(
            date, "(.+)/(.+)/(.+)", "20\\3-\\1-\\2"
        ) %>%
            ymd()
    ) %>%
    dplyr::rename(
        Country = `Country/Region`
    ) %>%
    dplyr::mutate(
        Country = case_when(
            `Province/State` == "Hubei" ~ "China (Hubei)",
            `Province/State` == "Hong Kong" ~ "Hong Kong",
            grepl("China", Country) & !`Province/State` %in% c("Hubei", "Hong Kong") ~ "China (Other)",
            Country == "Korea, South" ~ "South Korea",
            !grepl("China", Country) ~ Country
        )
    ) %>%
    dplyr::filter(
        !is.na(confirmed)
    ) %>%
  dplyr::select(-Lat, -Long)

govTable <- tibble()
govUrl <- "https://www.health.gov.au/news/health-alerts/novel-coronavirus-2019-ncov-health-alert/coronavirus-covid-19-current-situation-and-case-numbers"
if (url.exists(govUrl)) {
  govDate <- govUrl %>%
    read_html() %>%
    html_nodes("body") %>%
    xml_find_all("//p[contains(@class, 'au-callout')]") %>%
    html_text() %>%
    .[[1]] %>% 
    str_split("\n") %>%
    .[[1]] %>% 
    str_trim() %>% 
    .[. != ""] %>% 
    str_replace_all(".+ ([0-9]+) (March|April) (2020).+", "\\3-\\2-\\1") %>% 
    str_replace_all("March", "03") %>%
    str_replace_all("April", "04") %>%
    ymd()
  govLines <- govUrl %>%
    read_html() %>%
    html_nodes("body") %>%
    xml_find_all("//div[contains(@class, 'health-table')]") %>% 
    html_text() 
  govTable <- govLines %>% 
    str_split("\n") %>% 
    .[[1]] %>% 
    str_trim() %>% 
    .[. != ""] %>% 
    matrix(ncol = 2, byrow = TRUE) %>%
    set_colnames(.[1,]) %>% 
    .[-1,] %>% 
    as_tibble() %>%
    dplyr::filter(
      Location %in% c(
        "Queensland",
        "Australian Capital Territory",
        "New South Wales",
        "Victoria",
        "Tasmania",
        "South Australia",
        "Western Australia",
        "Northern Territory"
      )
    ) %>% 
    rename(
      `Province/State` = Location, 
      confirmed = `Confirmed cases*`
    ) %>%
    mutate(
      confirmed = str_remove_all(confirmed, ",") %>% as.numeric(),
      date = govDate,
      Country = "Australia"
    )
}

guardianData <- fromJSON("https://interactive.guim.co.uk/docsdata/1q5gdePANXci8enuiS4oHUJxcxC13d6bjMRSicakychE.json")$sheets 
guardTable <- guardianData$`latest totals` %>%
  as_tibble() %>%
  dplyr::select(
    `Province/State` = `Long name`,
    confirmed = `Confirmed cases (cumulative)`,
    date = `Last updated`,
    deaths = Deaths
  ) %>%
  mutate(
    Country = "Australia",
    date = ymd(date),
    confirmed = as.numeric(confirmed),
    deaths = as.numeric(deaths),
    deaths = ifelse(is.na(deaths), 0, deaths)
  ) %>%
  dplyr::filter(`Province/State` != "National")

nswTests <- NA_real_
nswUrl <- "https://www.health.nsw.gov.au/Infectious/diseases/Pages/covid-19-latest.aspx"
nswTests <- nswUrl %>%
  read_html() %>%
  html_nodes("body")  %>%
  xml_find_all("//td[contains(@class, 'moh-rteTableFooterOddCol-6')]") %>%
  html_text() %>%
  .[[1]] %>%
  str_remove(",") %>%
  as.numeric()

qldTests <- NA_real_
qldUrl <- glue(
  "https://www.qld.gov.au/health/conditions/health-alerts/coronavirus-covid-19/current-status/current-status-and-contact-tracing-alerts#symptoms-testing"
)
if (url.exists(qldUrl)){
  qldTests <- qldUrl %>%
    read_html() %>%
    html_nodes("body") %>%
    xml_find_all("//table[contains(@id, 'table59454')]") %>%
    html_text() %>%
    str_replace_all(
      ".+Samples testedTotal([0-9,]+)Tests.+",
      "\\1"
    ) %>%
    str_remove_all(",") %>%
    as.numeric()
}

vicRecov <- vicTests <- NA_real_
vicUrl <- glue(
  "https://www.dhhs.vic.gov.au/coronavirus-update-victoria-{str_trim(format(dt, '%e-%B-%Y'))}"
) %>%
  str_to_lower()
if (url.exists(vicUrl)){
  vicTxt <- vicUrl %>%
    read_html() %>%
    html_nodes("body") %>%
    xml_find_all("//div[contains(@class, 'field field--name-field-dhhs')]") %>%
    html_text() %>%
    str_split("\n") %>%
    .[[1]] 
  vicTests <- vicTxt %>%
    str_subset("test(ed|s)") %>%
    str_replace_all(
      ".+than ([0-9,]+) .*test(ed|s).+",
      "\\1"
    ) %>%
    str_remove_all(",") %>%
    as.numeric()
  vicRecov <- vicTxt %>%
    str_subset("recovered") %>%
    str_replace_all(".+ ([0-9,]+) people have recovered.", "\\1") %>%
    str_remove(",") %>%
    as.numeric()
}

waTests <- NA_real_
waRecov <- NA_real_
# waUrl <- glue(
#   "https://ww2.health.wa.gov.au/Media-releases/2020/COVID19-update-{str_replace(format(dt, '%d-%B-%Y'), '^0', '')}"
# )
waUrl <- "https://ww2.health.wa.gov.au/Articles/A_E/Coronavirus/COVID19-statistics"
waTxt <- waUrl %>%
  read_html() %>%
  xml_find_all("//table") %>%
  .[[1]] %>%
  html_text() %>%
  str_split("\n") %>%
  .[[1]] %>%
  str_trim() %>%
  .[. != ""] %>%
  str_remove(",") %>%
  matrix(ncol = 2, byrow = TRUE) %>%
  set_colnames(.[1,]) %>%
  .[-1,] %>%
  as_tibble() %>%
  mutate(Total = as.numeric(Total))
waTests <- waTxt %>%
  dplyr::filter(str_detect(Cases, "ive")) %>%
  summarise(n = sum(Total)) %>%
  .[["n"]]
waRecov <- dplyr::filter(waTxt, Cases == "Recovered")$Total

saRecov <- saTests <- NA_real_
saUrl <- "https://www.covid-19.sa.gov.au/home/dashboard"
saTxt <- saUrl  %>%
  read_html() %>%
  xml_find_all("//div[contains(@class, 'accordion__content')]") %>%
  .[[1]] %>%
  html_text() %>%
  str_split("\r\n") %>%
  .[[1]] %>%
  str_trim() %>%
  .[. != ""]
saTests <- saTxt %>%
  str_replace_all(".+have been.* ([0-9, ]+) COVID-19 laboratory tests.+", "\\1") %>%
  str_remove_all("[, ]") %>%
  as.numeric()
saRecov <- saTxt %>%
  str_replace(".+\\..+ ([0-9]+) .*recovered.+", "\\1") %>%
  as.numeric()

actRecov <- actTests <- NA_real_
actUrl <- "https://www.health.act.gov.au/about-our-health-system/novel-coronavirus-covid-19"
if (url.exists(actUrl)){
  actTxt <- actUrl %>%
    read_html() %>%
    html_nodes("body") %>%
    xml_find_all(
       "//div[contains(@class, 'spf-article-card--tabular')]"
    ) %>% 
    html_text() %>% 
    str_split("\r\n") %>% 
    .[[4]] %>% 
    str_trim() %>% 
    setdiff(y = "") %>% 
    str_remove_all("[,*]") %>%
    matrix(ncol = 2, byrow = TRUE) 
  actTests <- actTxt %>% 
    .[!str_detect(.[,1], "recovered"), 2] %>%
    as.numeric() %>% 
    sum()
  actRecov <- actTxt %>%
    .[str_detect(.[,1], "recov"), 2] %>%
    as.numeric()
}

tasTests <- NA_real_
tasUrl <- glue(
  "https://www.dhhs.tas.gov.au/news/2020/coronavirus_update_{format(dt, '%d_%B_%Y')}"
) %>%
  str_to_lower()
if (url.exists(tasUrl)) {
  suppressWarnings(
    tasTests <- tasUrl %>%
      read_html() %>%
      html_nodes("body") %>%
      xml_find_all(
        "//div[contains(@id, 'content_div_117564')]"
      ) %>%
      html_text() %>%
      str_split("\r\n") %>%
      .[[1]] %>%
      str_trim() %>%
      .[. != ""] %>%
      str_replace_all(".+Tasmania has conducted ([0-9]+) tests.+", "\\1") %>%
      as.numeric()
  )
} 

ntRecov <- ntTests <- NA_real_
ntUrl <- "https://coronavirus.nt.gov.au/"
ntTxt <- ntUrl %>%
  read_html() %>%
  html_nodes("body") %>%
  xml_find_all(
    "//p[contains(@class, 'mt-2')]"
  ) %>%
  html_text() %>%
  .[1] 
ntTests <- ntTxt %>%
  str_replace_all("^([0-9,]+) (people|test).+", "\\1") %>%
  str_remove_all(",") %>%
  as.numeric()
ntRecov <- ntTxt %>%
  str_replace(".+conducted([0-9,]+) .+recovered", "\\1") %>%
  as.numeric()

auRec <-c(
  "Australian Capital Territory" = actRecov,
  "Northern Territory" = ntRecov,
  "South Australia" = saRecov,
  "Western Australia" = waRecov,
  "Victoria" = vicRecov
) %>%
  enframe("Province/State", "recovered") %>%
  mutate(
    date = ymd(dt),
    Country = "Australia"
  ) %>%
  bind_rows(
    guardianData$updates %>% 
      as_tibble() %>%
      dplyr::select(State, Date, recovered = `Recovered (cumulative)`) %>%
      mutate(
        date = str_replace_all(Date, "([0-9]+)/([0-9]+)/2020", "2020-\\2-\\1") %>% ymd(),
        recovered = as.numeric(recovered),
        `Province/State` = case_when(
          State == "ACT" ~ "Australian Capital Territory" ,
          State == "NSW" ~ "New South Wales",
          State == "QLD" ~ "Queensland",
          State == "VIC" ~ "Victoria",
          State == "SA" ~ "South Australia",
          State == "WA" ~ "Western Australia",
          State == "TAS" ~ "Tasmania",
          State == "NT" ~ "Northern Territory"
        ),
        Country = "Australia"
      ) %>% 
      dplyr::filter(!is.na(recovered)) 
  ) %>%
  bind_rows(
    here::here("recovered") %>%
      list.files(full.names = TRUE, pattern = "recovered.+tsv", ) %>%
      sort() %>%
      .[[length(.)]] %>%
      read_tsv()
  ) %>%
  group_by(`Province/State`, Country) %>%
  mutate(
    recovered = max(recovered)
  ) %>%
  ungroup() %>%
  arrange(desc(recovered), desc(date)) %>%
  distinct(`Province/State`, .keep_all = TRUE) %>%
  mutate(date = max(date)) %>%
  dplyr::select(`Province/State`, Country, date, recovered)
if (nrow(auRec) > 0 ){
  auRec %>%
    write_tsv(
      here::here(glue("recovered/recovered_{dt}.tsv"))
    )
}

prevDateAU <- confirmed %>%
  dplyr::filter(Country == "Australia") %>%
  summarise_at(vars(date), max) %>%
  .[["date"]]
latestAU <- fromJSON("https://api.infotorch.org/api/covid19/totals/") %>%
  as_tibble() %>%
  rename(
    `Province/State` = state_long,
    date = date_updated
  ) %>%
  mutate(
    date = ymd(date),
    Country = "Australia",
    `Province/State` = str_replace_all(
      `Province/State`, "ACT", "Australian Capital Territory"
    )
  ) %>%
  dplyr::filter(
    !str_detect(`Province/State`, "Total"),
    `Province/State` != "Australia"
  ) %>%
  dplyr::select(
    `Province/State`, Country, date, 
    confirmed, deaths
  ) %>%
  bind_rows(govTable) %>% 
  bind_rows(guardTable) %>%
  # dplyr::filter(date > prevDateAU) %>%
  group_by(`Province/State`) %>% 
  summarise(
    Country = unique(Country),
    date = max(date),
    confirmed = max(confirmed), 
    deaths = max(deaths, na.rm = TRUE)
  ) %>% 
  mutate(
    deaths = ifelse(is.finite(deaths), deaths, NA)
  ) %>%
  left_join(auRec)
## The API & Govt website seems not to keep on top of NT & TAS cases so well
## Future versions of this script may need to manually update these via webscraping
## Here, just ensure that we don't have a decrease
oldTas <- max(dplyr::filter(confirmed, `Province/State` == "Tasmania")$confirmed)
newTas <- max(dplyr::filter(latestAU, `Province/State` == "Tasmania")$confirmed)
oldNT <- max(dplyr::filter(confirmed, `Province/State` == "Northern Territory")$confirmed)
newNT <- max(dplyr::filter(latestAU, `Province/State` == "Northern Territory")$confirmed)
if (newTas < oldTas | newNT < oldNT) {
  latestAU %<>%
    mutate(
      confirmed = case_when(
        `Province/State` == "Tasmania" ~ max(newTas, oldTas),
        `Province/State` == "Northern Territory" ~ max(newNT, oldNT),
        !`Province/State` %in% c("Tasmania", "Northern Territory") ~ confirmed
      )
    )
}
# Update the AU records if there are any timepoints
# on days beyond the last value provided by JHU
updateAU <- max(latestAU$date) > prevDateAU
if (updateAU) {
  confirmed %<>%
    bind_rows(
      dplyr::select(latestAU, any_of(colnames(.)))
    )
  n <- confirmed %>%
    dplyr::filter(Country == "Australia", date == prevDateAU) %>%
    bind_rows(latestAU) %>% 
    group_by(`Province/State`) %>% 
    summarise(n = length(unique(confirmed))) %>% 
    dplyr::filter(n > 1) %>%
    nrow()
}

recovered <- url("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv") %>%
  read_csv() %>%
  pivot_longer(
    cols = ends_with("20"),
    names_to = "date",
    values_to = "recovered"
  )  %>%
  mutate(
    date = str_replace_all(
      date, "(.+)/(.+)/(.+)", "20\\3-\\1-\\2"
    ) %>%
      ymd()
  ) %>%
  dplyr::rename(
    Country = `Country/Region`
  ) %>%
  dplyr::mutate(
    Country = case_when(
      `Province/State` == "Hubei" ~ "China (Hubei)",
      `Province/State` == "Hong Kong" ~ "Hong Kong",
      grepl("China", Country) & !`Province/State` %in% c("Hubei", "Hong Kong") ~ "China (Other)",
      grepl("Korea, South", Country) ~ "South Korea",
      !grepl("China", Country) ~ Country
    )
  ) %>%
  dplyr::select(-Lat, -Long)
if (updateAU){
  recovered %<>%
    bind_rows(
      latestAU %>%
        dplyr::select(any_of(colnames(recovered)))
    ) %>%
    arrange(
      `Province/State`, Country, date
    ) %>%
    mutate(recovered = zoo::na.locf(recovered)) %>%
    group_by(
      `Province/State`, Country, date
    ) %>%
    summarise_at(vars(recovered), max)
}

deaths <- url("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv") %>%
  read_csv() %>%
  pivot_longer(
    cols = ends_with("20"),
    names_to = "date",
    values_to = "deaths"
  )  %>%
  mutate(
    date = str_replace_all(
      date, "(.+)/(.+)/(.+)", "20\\3-\\1-\\2"
    ) %>%
      ymd()
  ) %>%
  dplyr::rename(
    Country = `Country/Region`
  ) %>%
  dplyr::mutate(
    Country = case_when(
      `Province/State` == "Hubei" ~ "China (Hubei)",
      `Province/State` == "Hong Kong" ~ "Hong Kong",
      grepl("China", Country) & !`Province/State` %in% c("Hubei", "Hong Kong") ~ "China (Other)",
      Country == "Korea, South" ~ "South Korea",
      !grepl("China", Country) ~ Country
    )
  ) %>%
  dplyr::select(-Lat, -Long)
if (updateAU){
  deaths %<>%
    bind_rows(
      dplyr::select(latestAU, any_of(colnames(.)))
    ) %>%
    mutate(deaths = zoo::na.locf(deaths))
}

wikiPops <- read_tsv("wikiPops.tsv")

Data for confirmed cases, recoveries and fatalities was primarily sourced from Johns Hopkins University (https://coronavirus.jhu.edu/), using the datasets provided at https://github.com/CSSEGISandData/COVID-19. JHU data is updated every 24 hours at approximately 23:59(UTC), which is about 9:30AM in Adelaide.

Live hourly updates for Australia are available from https://covid-19-au.github.io/ for those who would like an up to the minute breakdown of confirmed cases. Numbers used for generation of this page are updated periodically throughout the day using values provided by https://api.infotorch.org/api/covid19/totals/, The Guardian and those at www.health.gov.au. The Australian numbers used below do contain confirmed cases reported since the JHU dataset was updated. As of 12:28, 4 of 8 states/territories have released updated COVID-19 figures since the overnight JHU release. The current Australia total is 5,668 using these data sources, and any corrections to these will be made as soon as is possible.

Population sizes were obtained from 2019 UN Estimates. Given the disparity of infection within China, China was broken into Hubei Province and the rest of China, with Hong Kong and Taiwan already being considered as separate countries in all datasets. Population estimates for Hubei Province were taken from the 2018 estimates given by Statista.com and this is likely to be a very slight underestimate.

Confirmed cases of COVID-19 as provided by the Chinese Government have been discussed elsewhere as unusual, and data appears potentially unreliable. In this analysis, discussions regarding accurate Chinese reporting are not considered further and data is simply presented as supplied by JHU.

However, all countries are likely to contain many unreported cases given the incomplete testing regimes in place for most countries. Similarly, reporting in many countries may have features that cause concerns regarding data integrity and this makes comparison across countries difficult. Information on recovered cases is also difficult to accurately obtain due inconsistent methods for considering a case as recovered, and lack of reporting for these cases in many jurisdictions.

International Data

startingPoint <- 4
minPop <- 4e6

For this section, most data is presented relative to population size. Growth in infection rates is only shown after the point at which the cumulative confirmed infection rate breached 4 confirmed cases / million. This equates to about 101 confirmed cases within Australia, and is broadly comparable to the “Days since passing 100 confirmed cases” commonly shown elsewhere.

The exception to this is the Daily Increase Vs Confirmed Cases plot which uses raw case numbers and 100 confirmed cases as the point for beginning the plot.

Table of Most Impacted Countries

Confirmed cases in this table are effectively the cumulative, confirmed incidence rate. Recovered patients and those who have passed away are still included in these numbers. At the time of report preparation, the total number of confirmed cases, including all countries for which data has been made available, is 1,197,523.

confirmed %>%
  group_by(Country, date) %>%
  summarise(confirmed = sum(confirmed)) %>%
  ungroup() %>%
  group_by(Country) %>%
  dplyr::filter(
    date == max(date),
  ) %>%
  ungroup() %>%
  inner_join(wikiPops) %>%
  mutate(
    rate = 1e6*confirmed / Population,
    occurrence = Population / confirmed,
    Population = round(Population, -3) / 1e6
  ) %>%
  dplyr::filter(rate >= 1) %>%
  arrange(desc(rate)) %>%
  rename_at(vars(everything()), str_to_title) %>%
  dplyr::select(
    Continent, Region, Country,
    Date, Confirmed, 
    `Population (millions)` = Population, 
    `Rate (Cases per Million)` = Rate,
    Occurrence
  ) %>%
  datatable(
    options = list(
      pageLength = 25, 
      autoWidth = TRUE,
      searchCols = list(
        NULL, NULL, NULL, NULL, NULL,
        list(
          search = glue(
            '{minPop/1e6 + 0.001} ... {max(wikiPops$Population/1e6)}'
          )
        ),
        NULL
      )
    ),
    filter = 'top',
    class = "stripe",
    rownames = FALSE,
    caption = paste(
      "The most impacted countries studied here and shown as a proportion of total population.",
      "The initial filter is set so that only countries with a population greater than", comma(minPop), "are shown.",
      "All fields are searchable and sortable.",
      "To filter numeric columns, either use the slider or enter the values in the form 'min ... max'.",
      "To filter text columns, partial matching used in a case-insensitive manner.",
      "Populations have been rounded to the nearest thousand to make reading values easier.",
      "'Rate' represents the latest confirmed infection rate as cumulative cases per million people, whilst 'Occurrence' represents the number of people expected before one case is found, assuming an even distribution amongst the population.",
      "In other words, one in every 'Occurrence' people within the population have been confirmed to have contracted COVID-19.",
      "Occurrence is inversely proportional to Rate.",
      "No adjustment has been made in this table for patients who have recovered or passed away.",
      "Whilst the virus spreads with no regard to population size, the rate as shown here indicates the degree of stress which each country's health-care system is likely to be experiencing.",
      "Several countries shown here have not attracted much media attention due lower case numbers than China and Italy, but are likely to be experiencing significant duress.",
      "Continent and Region information is as provided by the UN classifications."
    )
  ) %>%
  formatCurrency(
    columns = c("Population (millions)"),
    currency = "", 
    digits = 3,
    mark = ","
  ) %>%
  formatCurrency(
    columns = c("Confirmed", "Rate (Cases per Million)", "Occurrence"),
    currency = "", 
    digits = 0,
    mark = ","
  )

Daily Increase Vs Confirmed Cases

An alternate viewpoint on the data is to remove time and inspect the relationship between the daily increase in cases and the total number of cases. When this relationship ceases it’s near linear relationship, this can be a sign the control measures have begun to take effect.

sma <- function(x, n = 7){
    
    if (length(x) < n) return(rep(NA, length(x)))
    ind <- vapply(seq(1, length(x) - (n - 1)), seq, numeric(n), length.out = n)
    c(
        rep(NA, n - 1),
        colMeans(matrix(x[ind], nrow = n), na.rm = TRUE)
    )
}

minPop <- 10e6
minCases <- 100
minDays <- 15
plotIncVConf <- confirmed %>%
  group_by(Country, date) %>%
  summarise_at(vars(confirmed), sum) %>%
  dplyr::filter(confirmed > 0) %>%
  split(f = .$Country) %>%
  lapply(function(x, n = 7){
    x %>%
      mutate(
        d = c(0, diff(confirmed))
      ) %>%
      mutate_at(
        vars(confirmed, d), sma, n = n
      )
  }
  ) %>%
  bind_rows() %>%
  inner_join(wikiPops) %>%
  dplyr::filter(confirmed > minCases) %>%
  mutate_at(vars(confirmed, d), round, 2) %>%
  mutate(n = n()) %>%
  dplyr::filter(
    n > minDays,
    Population > minPop,
    d > 0
  ) %>%
  rename_all(str_to_title) %>%
  rename(`Average Daily Increase` = D) %>%
  mutate(Population = comma(round(Population, -3))) %>%
  ggplot(aes(Confirmed, `Average Daily Increase`, colour = Country, label = Date, key = Population)) +
  geom_hline(
    aes(yintercept = `Average Daily Increase`),
    data = . %>%
      dplyr::filter(Country == "Australia", Date == max(Date)),
    linetype = 3,
    colour = "blue",
    size = 1/3
  ) +
  geom_vline(
    aes(xintercept = Confirmed),
    data = . %>%
      dplyr::filter(Country == "Australia", Date == max(Date)),
    linetype = 3,
    colour = "blue",
    size = 1/3
  ) +
  geom_point() +
  geom_line() +
  scale_x_log10(labels = comma) +
  scale_y_log10(labels = comma) +
  labs(
    x = "Confirmed Cases",
    y = 'Average Daily Increase'
  )

capIncVConf <- glue(
  "*Daily Increase in Cases plotted against Confirmed Cases.
  These two values are directly proportional until interventions are successful, at which point the proportional relationship changes, as evidenced by a sudden downwards turn.
  Values shown are 7-day moving averages in order to minimise the impact of day-to-day variation.
  Countries are only shown from the point at which the moving average exceeds {minCases} cases, and have exceeded this value for > {minDays} days.
  Data is additionally restricted to countries with a population > {comma(minPop)}, for greater clarity.
  Australia's position is indicated by the intersection of the dashed blue lines.*
  "
)

ggplotly(plotIncVConf)

Daily Increase in Cases plotted against Confirmed Cases. These two values are directly proportional until interventions are successful, at which point the proportional relationship changes, as evidenced by a sudden downwards turn. Values shown are 7-day moving averages in order to minimise the impact of day-to-day variation. Countries are only shown from the point at which the moving average exceeds 100 cases, and have exceeded this value for > 15 days. Data is additionally restricted to countries with a population > 10,000,000, for greater clarity. Australia’s position is indicated by the intersection of the dashed blue lines.

Cumulative Incidence Rates

ausDays <- confirmed %>% 
  dplyr::filter(Country == "Australia") %>%
  group_by(Country, date) %>%
  summarise_at(vars(confirmed), sum) %>%
  left_join(wikiPops) %>%
  mutate(rate = 1e6 * confirmed / Population) %>%
  dplyr::filter(rate > startingPoint) %>%
  nrow()
minDays <- ausDays - 5
# Use Singapore as that has the longest dataset besides Hubei
nDays <- confirmed %>%
  dplyr::filter(Country == "Singapore") %>% 
  group_by(Country, date) %>%
  summarise(confirmed = sum(confirmed)) %>%
  ungroup() %>%
  left_join(wikiPops) %>%
  mutate(
    rate = 1e6*confirmed / Population
  ) %>% 
  dplyr::filter(rate > startingPoint) %>% 
  nrow() %>%
  subtract(1)

refRate <- c(2, 4, 8)
minPop <- 5e6
p <- confirmed %>%
  group_by(Country, date) %>%
  summarise(confirmed = sum(confirmed)) %>%
  ungroup() %>%
  inner_join(
    dplyr::filter(wikiPops, Population > minPop)
  ) %>%
  mutate(
    rate = 1e6*confirmed / Population
  ) %>%
  dplyr::filter(
    rate > startingPoint
  ) %>%
  group_by(Country) %>%
  mutate(
    days = date - min(date)
  ) %>%
  dplyr::filter(
    max(days) >= minDays
  ) %>%
  ungroup() %>%
  mutate(
    days = as.integer(days),
    rate = round(rate, 2)
  ) %>%
  dplyr::filter(days <= nDays) %>%
  arrange(date) %>%
  mutate(Country = fct_inorder(Country)) %>%
  rename_all(str_to_title) %>%
  ggplot(
    aes(Days, Rate, colour = Country, Date = Date, Confirmed = Confirmed)
  ) +
  geom_segment(
    aes(x, y, xend = xmax, yend = ymax),
    data = tibble(
      x = 0,
      y = startingPoint,
      xmax = c(16, 32, nDays ),
      ymax = startingPoint*2^(xmax / refRate)
    ),
    inherit.aes = FALSE,
    colour = "grey70",
    linetype = 2
  ) +
  geom_text(
    aes(xmax, ymax, label = label),
    data = tibble(
      xmax = c(16, 32, nDays - 2),
      ymax = startingPoint*2^(xmax / refRate),
      label = glue("Doubling in\n {refRate} days")
    ),
    colour = "grey70",
    inherit.aes = FALSE
  ) +
  geom_vline(
    aes(xintercept = Days),
    data = . %>%
      dplyr::filter(Country == "Australia") %>%
      dplyr::filter(Days == max(Days)),
    linetype = 3,
    colour = "blue",
    size = 1/3
  ) +
  geom_hline(
    aes(yintercept = Rate),
    data = . %>%
      dplyr::filter(Country == "Australia") %>%
      dplyr::filter(Days == max(Days)),
    linetype = 3,
    colour = "blue",
    size = 1/3
  ) +
  geom_point() +
  geom_line() +
  scale_x_continuous(
    expand = expand_scale(mult = c(0, 0.05)),
  ) +
  scale_y_log10(
    expand = expand_scale(mult = c(0, 0.05))
  ) +
  xlab(
    paste(
      "Days since passing", 
      startingPoint,
      "confirmed cases/million"
    )
  ) +
  ylab("Confirmed Cumulative Infection Rate (cases/million)")

ggplotly(
  p, 
  tooltip = c(
    "Days", "Rate", "Country", "Date", "Confirmed"
  )) 

COVID-19 Confirmed Cumulative Infection Rate for countries which have exceeded 4 confirmed cases/million for 22 or more days, and with populations greater than 5,000,000. Data is only shown for the first 60 calendar days since passing 4 confirmed cases/million. Note that from the day records begin in this dataset (2020-01-22), the confirmed infection rate in Hubei was 7.5 confirmed cases/million. The blue dashed line indicates Australia’s current position on this figure. Diagonal grey lines indicate a doubling in the infection rate every 2, 4, or 8 days. To hide a country, click on the country in the plot legend. Clicking again on the country in the legend will restore the data within the plot. Countries are shown in order of the date at which they passed the 4 confirmed case/million mark. Due to the large number of countries shown, you may need to scroll through the legend. Regions of the plot are also able to be zoomed interactively. Please note the y-axis is shown on the logarithmic scale, so that a series of points which appear to be diagonal will indicate exponential growth. The flatter the line, the slower the growth and a perfectly horizontal line would indicate zero growth, or no new confirmed cases.

All figures and tables presented here simply aim to show an alternative viewpoint on the data. Every way to view COVID-19 data will mask important features, and the values shown here do not take into account vital factors such as population density, variability of infection across regions within countries, social culture and demographics. Many countries may not be directly comparable for a combination of the above factors. It is simply to view the data through the lens of a country’s population size using a value which should be easily interpretable.

In the above plot:

• Only countries are shown where the infection rate has surpassed 4 cases/million for more than 22 days. This choice seemed to be a useful way of enabling the comparison of COVID-19 progression across countries after scaling for population size
• Only South Korea,and China appear to have brought the virus under some level of control. The infection rate currently sits around 198.3 per million for South Korea and 1145.9 per million for Hubei Province. This can also be thought of as 0.0198% and 0.1146% of the total populations respectively. Note that despite these apparently low percentages, the actual numbers of deaths and hospitalisations are extremely non-trivial and the real impacts of these numbers cannot be treated lightly.
• Japan’s early measures also appear to be highly effective, however the recent growth in infections within Singapore and Hong Kong should serve as a caution
• US numbers may also be unreliable given the low rates of testing and numerous anecdotal stories of suspected infections

Currently Active Infection Rates

As an alternative viewpoint, the numbers of recovered and deceased patients have been removed from the numbers of confirmed infections plot to provide an estimate of the currently active infections. Importantly, information regarding recovered cases is likely to be the least reliable of reported values as many regions do not report updated numbers for several consecutive days. Additionally many regions do not report recovered cases as the criteria for considering a person to have recovered as currently unclear.

p2 <- confirmed %>%
  dplyr::filter(
    confirmed > 0,
    Country != "China (Other)"
  ) %>%
  left_join(deaths) %>%
  group_by(Country, date) %>%
  summarise_at(
    vars(confirmed, deaths), sum
  ) %>%
  left_join(
    recovered %>%
      group_by(Country, date) %>%
      summarise_at(vars(recovered), sum)
  ) %>%
  mutate(
    active = confirmed - deaths - recovered
  ) %>%
  dplyr::filter(!is.na(active)) %>%
  inner_join(
    wikiPops %>%
      dplyr::filter(Population > minPop)
  ) %>%
  mutate(rate = 1e6 * active / Population) %>%
  dplyr::filter(rate > startingPoint) %>%
    group_by(Country) %>%
  mutate(
    days = date - min(date)
  ) %>%
  dplyr::filter(max(days) > minDays) %>%
  ungroup() %>%
  mutate(
    days = as.integer(days),
    rate = round(rate, 2)
  ) %>%
  dplyr::filter(days <= nDays) %>%
  arrange(date) %>%
  mutate(Country = fct_inorder(Country)) %>%
  rename_all(str_to_title) %>%
  ggplot(
    aes(
      x = Days, y = Rate, colour = Country, 
      Date = Date, Active = Active,
      Confirmed = Confirmed, Recovered = Recovered,
      Deaths = Deaths
      )
  ) +
  geom_segment(
    aes(x, y, xend = xmax, yend = ymax),
    data = tibble(
      x = 0,
      y = startingPoint,
      xmax = c(21, 42, nDays),
      ymax = startingPoint + 2^(xmax/ c(2, 4, 8))
    ),
    inherit.aes = FALSE,
    colour = "grey70",
    linetype = 2
  ) + 
  geom_text(
        aes(xmax, ymax, label = label),
    data = tibble(
      xmax = c(21, 42, nDays),
      rate = c(2, 4, 8),
      ymax = 2^(xmax/ rate),
      label = glue("Double in\n {rate} days")
    ),
    colour = "grey70",
    inherit.aes = FALSE
  ) +
  geom_vline(
    aes(xintercept = Days),
    data = . %>%
      dplyr::filter(Country == "Australia") %>%
      dplyr::filter(Days == max(Days)),
    linetype = 3,
    colour = "blue",
    size = 1/3
  ) +
  geom_hline(
    aes(yintercept = Rate),
    data = . %>%
      dplyr::filter(Country == "Australia") %>%
      dplyr::filter(Days == max(Days)),
    linetype = 3,
    colour = "blue",
    size = 1/3
  ) +
  geom_point() +
  geom_line() +
  scale_x_continuous(
    expand = expand_scale(mult = c(0, 0.05)),
  ) +
  scale_y_log10(
    expand = expand_scale(mult = c(0, 0.05))
  ) +
  xlab(
    paste(
      "Days since passing", 
      startingPoint, 
      "confirmed active cases/million"
    )
  ) +
  ylab("Confirmed Active Infection Rate (cases/million)")

ggplotly(p2) 

Confirmed active cases of COVID-19 for countries where the confirmed infection rate has exceeded 4 confirmed active cases/million for more than 60 calendar days. Only countries with a population greater than 5,000,000 are shown for better visualisation. Due to difficulties introduced by the currently reported low active infection rate outside Hubei province, data from China has been excluded from this plot, with the exception of Hubei and Hong Kong. The blue dashed lines indicates Australia’s current position on this figure. To hide a country, click on the country in the plot legend. Clicking again on the country in the legend will restore the data within the plot. Countries are shown in order of the date at which they passed the 4 confirmed active case/million mark. Due to the number of countries shown, you may need to scroll through the legend. Regions of the plot are also able to be zoomed interactively. Please note the y-axis is shown on the logarithmic scale, so that a series of points which appear to be diagonal will indicate exponential growth/decay. Given the different starting point to the previous plot, data will generally be shown for fewer time-points.

Notable features of this perspective are:

• Active infection rates in South Korea are currently either declining or stable
• Despite strict control measures, active infection rates are exponentially increasing in both Hong Kong and Singapore

Dimensional Reduction

In order to summarise which countries are the most similar to each other, a Principal Component Analysis was performed. This enables the multi-dimensional data of the above plots to summarised in two dimensions. As missing data cannot be included in this analysis, several countries which are at earlier comparative time-points than Australia were omitted from this analysis.

nDays <- confirmed %>% 
  dplyr::filter(Country == "Australia") %>% 
  group_by(Country, date) %>%
  summarise(confirmed = sum(confirmed)) %>%
  ungroup() %>%
  inner_join(wikiPops) %>%
  mutate(
    rate = 1e6*confirmed / Population
  ) %>% 
  dplyr::filter(rate > startingPoint) %>% 
  nrow() %>%
  subtract(1)

minPop <- 2e6
pca <- confirmed %>%
  group_by(Country, date) %>%
  summarise(confirmed = sum(confirmed)) %>%
  ungroup() %>%
  inner_join(wikiPops) %>%
  mutate(
    rate = 1e6*confirmed / Population
  ) %>%
  dplyr::filter(
    rate > startingPoint,
    Population > minPop
  ) %>%
  group_by(Country) %>%
  mutate(
    days = date - min(date)
  ) %>%
  dplyr::filter(
    max(days) >= nDays
  ) %>%
  ungroup() %>%
  mutate(
    days = as.integer(days),
    rate = round(rate, 2)
  ) %>%
  dplyr::filter(days <= nDays) %>%
  dplyr::select(Country, rate, days) %>%
  pivot_wider(
    values_from = rate,
    names_from = days
  ) %>%
  as.data.frame() %>%
  column_to_rownames("Country") %>%
  as.matrix() %>%
  .[!rowAnyNAs(.),] %>%
  log() %>%
  prcomp() 

set.seed(101)
pca$x %>%
  as.data.frame() %>%
  rownames_to_column("Country") %>%
  mutate(
    Cluster = cbind(PC1, PC2) %>%
      apply(2, function(x){
        (x - mean(x)) / sd(x)
      }) %>%
      kmeans(centers = k, nstart = 10) %>%
      .[["cluster"]] %>%
      as.factor()
  ) %>%
  ggplot(aes(PC1, PC2, label = Country, colour = Cluster)) +
  geom_point() +
  geom_text_repel(
    show.legend = FALSE
  ) +
  stat_ellipse(
    aes(fill = Cluster),
    colour = NA,
    geom = "polygon",
    alpha = 0.1
    ) +
  xlab(
    paste0(
      "PC1 (", pcPercent[['PC1']],")"
      )
  ) +
    ylab(
    paste0(
      "PC2 (", pcPercent[['PC2']], ")"
      )
  ) +
  theme(
    legend.position = "none"
  )

Dimensional reduction showing which countries infection rates are the most similar to each other at the 26 day time point, after passing 4 confirmed cases/million. This may or may not be indicative of future spread within the population. The value 26 days was simply chosen as this marks how long since Australia passed this threshold. Countries which have not progressed beyond 4 confirmed cases/million for 26 days or more are not shown. Countries with populations < 2,000,000 are also excluded. Clusters were generated using k-means, manually specifying k = 4 and are not definitive. Principal Component 1, on the x-axis, corresponds to the greatest source of variability within the data (91.8%), and countries which appear near each other along this axis can be assumed to be showing similar growth in confirmed infection rates at this time point. Separation along the y-axis is less significant, but also worthy of note, as this represents 6.4% of variability within the data. At this early point, Australia’s cumulative, confirmed infection rate is diverging from the cluster of countries which have responded well and is becoming more similar to Israel, the UK and other poor responding countries. This is suggestive that the early measures instituted in Australia may have been inadequate.

Fatality, Recovery and Active Infection Rates

All rates presented in this section do not take population size into account, but simply look at the rates of recovery and fatalities within each country’s cohort.

Summary of All Rates

fr <- confirmed %>%
  inner_join(deaths) %>%
  group_by(Country) %>%
  dplyr::filter(date == max(date)) %>%
  ungroup() %>%
  summarise(fr = sum(deaths) / sum(confirmed)) %>%
  .[["fr"]]
rr <- confirmed %>%
  group_by(Country, date) %>%
  summarise_at(vars(confirmed), sum, na.rm = TRUE) %>%
  left_join(
    recovered %>% 
      group_by(Country, date) %>%
      summarise_at(vars(recovered), sum, na.rm = TRUE)
  ) %>%
  dplyr::filter(date == max(date)) %>%
  ungroup() %>%
  summarise(rr = sum(recovered) / sum(confirmed)) %>%
  .[["rr"]]

Summarising all available data from all countries:

• The current fatality rate from all confirmed cases is 5.4%. This may be a function of under-reporting of true cases and is very likely to be an overestimate.
• The current recovery rate from all confirmed cases is 20.6%. Given the level of under-reporting this may also be highly inaccurate.
• At this point, 74.0% of all confirmed infections are considered as ‘active’.

minPop <- 5e6
p4 <- confirmed %>%
  dplyr::filter(
    confirmed > 0
  ) %>%
  left_join(deaths) %>%
  group_by(Country, date) %>%
  summarise_at(vars(confirmed, deaths), sum) %>%
  left_join(
    recovered %>%
      group_by(Country, date) %>%
      summarise_at(vars(recovered), sum)
  ) %>%
  ungroup() %>%
  inner_join(
    wikiPops %>% dplyr::filter(Population > minPop)
  ) %>%
  mutate(rate = 1e6 * confirmed / Population) %>%
  dplyr::filter(rate > startingPoint) %>% 
  arrange(date)%>% 
  mutate(Country = fct_inorder(Country)) %>%
  group_by(Country) %>%
  dplyr::filter(date == max(date)) %>%
  mutate(
    active = confirmed - recovered - deaths
  ) %>%
  mutate(
    active = 100*active / confirmed,
    recovered = 100*recovered / confirmed,
    fatalities = 100*deaths / confirmed
  ) %>%
  dplyr::filter(active < 100) %>%
  pivot_longer(
    cols = c(active, recovered, fatalities),
    names_to = "Status",
    values_to = "Percentage"
  ) %>%
  mutate(
    Status = str_to_title(Status),
    Status = factor(
      Status, 
      levels = c("Active", "Recovered", "Fatalities")
    ),
    Percentage = round(Percentage, 2)
  ) %>%
  mutate(confirmed = comma(confirmed)) %>%
  rename(Confirmed = confirmed) %>%
  ggplot(
    aes(
      Country, Percentage,
      fill = Status, cases = Confirmed
    )
  ) +
  geom_col() +
  scale_fill_manual(
    values = c(
      Active = "blue",
      Recovered = "green",
      Fatalities = "red"
      )
  ) +
  scale_y_continuous(expand = expand_scale(0, 0)) +
  coord_flip() +
  labs(x = c()) +
  theme(
    legend.position = "none"
  )

ggplotly(p4)

Fatality, Recovery and Active Infection rates for countries which have exceeded 4 confirmed cases / million, and with a population size > 5,000,000. Countries are shown in order of the date they passed 4 confirmed cases / million

Fatality Rate Vs Time

minDays <- 12
minPop <- 5e6
df <- confirmed %>%
  inner_join(deaths) %>%
  group_by(Country, date) %>%
  summarise_at(
    vars(confirmed, deaths),
    sum
  ) %>%
  inner_join(
    wikiPops %>% dplyr::filter(Population > minPop)
  ) %>%
  ungroup() %>%
  mutate(
    `Infection Rate` = 1e6 * confirmed / Population
  ) %>%
  dplyr::filter(
    `Infection Rate` > startingPoint
    ) %>%
  group_by(Country) %>%
  mutate(
    Days = date - min(date),
    n = n()
  ) %>%
  dplyr::filter(
    n > minDays,
    max(deaths, na.rm = TRUE) > 0
  ) %>%
  mutate(
    Rate = deaths / confirmed,
    `Fatality Rate` = percent(Rate, accuracy = 0.1),
    minusT = date - max(date)
  ) %>%
  ungroup() 
plotFr <- mutate(
  df,
  Country = factor(
    Country, 
    levels = df %>%
      dplyr::select(Country, minusT, Rate) %>%
      pivot_wider(
        id_cols = Country, 
        names_from = minusT, 
        values_from = Rate
      ) %>% 
      as.data.frame() %>%
      column_to_rownames("Country") %>% 
      dist() %>% 
      hclust() %>% 
      as.dendrogram() %>% 
      labels()
  )
) %>% 
  rename_all(str_to_title) %>%
  ggplot(
    aes(
      x = Days, y = Country, fill = Rate,
      conf = Confirmed, 
      deaths = Deaths,
      date = Date,
      label = `Fatality Rate`
      )
  ) +
  geom_raster() +
  geom_vline(
    aes(xintercept = Days + 0.5),
    data = . %>%
      dplyr::filter(Country == "Australia") %>%
      dplyr::filter(Date == max(Date)),
    linetype = 3,
    size = 1/3,
    colour = "grey70"
  ) +
  scale_fill_viridis_c(
    option = "magma",
    breaks = seq(0, 0.12, by = 0.02)
  ) +
  scale_x_continuous(
    expand = expand_scale(0, 0),
    labels = abs
    ) +
  scale_y_discrete(expand = expand_scale(0, 0)) +
  labs(
    x = glue(
      "Days since passing {startingPoint} cases/million"
    ),
    y = c(),
    fill = "Fatality\nRate"
  ) +
  theme(
    panel.grid = element_blank()
  )

cpFr <- glue(
  "*Fatality Rate for confirmed cases after passing {startingPoint} confirmed cases / million.
  Only countries with {minDays} days of data beyond this time-point and a population size >{comma(minPop)} are shown.
  Country order is based on clustering using the most recent values.
  Countries with no recorded fatalities have been excluded for obvious reasons. 
  The dashed grey line indicates the time-point Australia is currently at.
  A clear trend of an increasing fatality rate with time is evident in many countries (e.g. Spain, France, Italy, UK), however, for some countries (e.g. Singapore) this rate appears relatively stable throughout the majority of days.
  The overall Fatality Rate for confirmed cases is currently ({percent(fr, accuracy = 0.1)}).*"
)

ggplotly(
  plotFr,
  tooltip = c(
    "Country", "Date", "Days", "Confirmed", "Deaths",
    "Fatality Rate"
    )
  )

Fatality Rate for confirmed cases after passing 4 confirmed cases / million. Only countries with 12 days of data beyond this time-point and a population size >5,000,000 are shown. Country order is based on clustering using the most recent values. Countries with no recorded fatalities have been excluded for obvious reasons. The dashed grey line indicates the time-point Australia is currently at. A clear trend of an increasing fatality rate with time is evident in many countries (e.g. Spain, France, Italy, UK), however, for some countries (e.g. Singapore) this rate appears relatively stable throughout the majority of days. The overall Fatality Rate for confirmed cases is currently (5.4%).

Recoveries Vs Fatalities Vs Time

minDays <- ausDays - 5
df <- confirmed %>%
  left_join(deaths) %>%
  group_by(Country, date) %>%
  summarise_at(vars(confirmed, deaths), sum) %>%
  left_join(
    recovered %>%
      group_by(Country, date) %>%
      summarise_at(vars(recovered), sum)
  ) %>%
  group_by(Country, date) %>%
  summarise_at(
    vars(confirmed, recovered, deaths),
    sum
  ) %>%
  ungroup() %>%
  right_join(
    wikiPops 
  ) %>%
  mutate(
    active = confirmed - deaths - recovered,
    `Infection Rate` = 1e6 * confirmed / Population
  ) %>%
  dplyr::filter(`Infection Rate` > startingPoint) %>%
  group_by(Country) %>%
  mutate(
    Days = date - min(date),
    ratio = round(log2(recovered / deaths), 2),
  ) %>%
  dplyr::filter(
    max(Days) > minDays
  ) %>%
  ungroup() %>%
  dplyr::filter(
    !is.infinite(ratio),
    !is.na(ratio)
  ) 
p6 <- df %>%
  mutate(
    Country = factor(
      Country, 
      levels = df %>%
        dplyr::select(Country, Days, ratio) %>%
        group_by(Country) %>% 
        mutate(minusT = as.integer(Days - max(Days))) %>% 
        pivot_wider(Country, minusT, values_from = ratio) %>% 
        as.data.frame() %>%
        column_to_rownames("Country") %>% 
        as.matrix %>% 
        .[matrixStats::rowCounts(.,value = NA) < ncol(.),] %>% 
        dist() %>% 
        hclust() %>%
        as.dendrogram() %>%
        labels()
    )
  ) %>%
  rename_all(str_to_title) %>%
  ggplot(
    aes(
      x = Days, y = Country, fill = Ratio,
      conf = Confirmed, 
      rec = Recovered,
      deaths = Deaths,
      date = Date,
    )
  ) +
  geom_raster() +
  scale_fill_gradient2(
    low = "red", high = "green", na.value = "grey80"
  ) +
  scale_x_continuous(
    expand = expand_scale(0, 0)
  ) +
  scale_y_discrete(
    expand = expand_scale(0, 0)
  ) +
  labs(
    x = glue(
      "Days since passing {startingPoint} cases/million"
    ),
    y = c(),
    fill = "Ratio"
  ) +
  theme(
    panel.grid = element_blank()
  )

cp6 <- glue(
  "*Comparison of recoveries and fatalities as a time course, beginning at the day cases exceeded {startingPoint} confirmed cases / million.
  Countries are only shown if data is available for more than {minDays} days since passing this threshold.
  The ratio of recoveries to fatalities is shown on the log~2~ scale, such that when fatalities outnumber recoveries, the colour become more red.
  Conversely if the recoveries outnumber the fatalities the colour will become more green.
  If these two numbers are approximately equal, the colour will be white.
  As the numbers are on the log~2~ scale, a ratio of 2, would indicate 2^2^ = 4 times the number of __recoveries to fatalities__.
  Similarly a ratio of -2 would indicate that fatalities outnumber recoveries by a factor of 2^2^ = 4.
  Ratios are only available once there is one or more confirmed recoveries in addition to one or more confirmed fatalities.*"
)

A time-course of the relationship between recoveries and fatalities is shown. Series of days which extend deeply into the red may possibly indicate where medical facilities are overloaded. However, it is possible that during these times more attention is focussed on keeping critically ill patients alive than confirming a recovery. It should be noted that many of the countries which are deep into the red at the current time are affluent Western countries. It is also worth noting that despite the high death toll in Italy, recoveries outnumber fatalities for the majority of this time-course.

ggplotly(p6)

Comparison of recoveries and fatalities as a time course, beginning at the day cases exceeded 4 confirmed cases / million. Countries are only shown if data is available for more than 22 days since passing this threshold. The ratio of recoveries to fatalities is shown on the log₂ scale, such that when fatalities outnumber recoveries, the colour become more red. Conversely if the recoveries outnumber the fatalities the colour will become more green. If these two numbers are approximately equal, the colour will be white. As the numbers are on the log₂ scale, a ratio of 2, would indicate 2² = 4 times the number of recoveries to fatalities. Similarly a ratio of -2 would indicate that fatalities outnumber recoveries by a factor of 2² = 4. Ratios are only available once there is one or more confirmed recoveries in addition to one or more confirmed fatalities.

Unscaled Fatality Rate

minDeaths <- 10
ausDays <- deaths %>%
  dplyr::filter(Country == "Australia") %>% 
  group_by(date) %>%
  summarise_at(vars(deaths), sum) %>% 
  dplyr::filter(deaths > 10) %>%
  nrow()
minDays <- ausDays - 2
minPop <- 5e6
deathPlot <-confirmed %>%
  left_join(deaths) %>%
  group_by(Country, date) %>%
  summarise_at(
    vars(confirmed, deaths),
    sum
  ) %>%
  ungroup() %>%
  dplyr::filter(
    deaths >= minDeaths,
    Country %in% dplyr::filter(wikiPops, Population > minPop)$Country
  ) %>%
  group_by(Country) %>%
  mutate(Days = as.integer(date - min(date))) %>%
  dplyr::filter(
    max(Days) > minDays
  ) %>%
  ungroup() %>%
  rename_all(str_to_title) %>%
  arrange(desc(Days)) %>%
  mutate(Country = fct_inorder(Country)) %>%
  ggplot(
    aes(
      x = Days, y = Deaths, colour = Country,
      label = Date, conf = Confirmed)
    ) +
    geom_vline(
    aes(xintercept = Days),
    data = . %>% 
      dplyr::filter(Country == "Australia") %>%
      slice(1),
    linetype = 3,
    colour = "blue"
  ) +
    geom_hline(
    aes(yintercept = Deaths),
    data = . %>% 
      dplyr::filter(Country == "Australia") %>%
      slice(1),
    linetype = 3,
    colour = "blue"
  ) +
  geom_point() +
  geom_line() +
  scale_y_log10() +
  labs(
    x = glue("Days Since Passing {minDeaths} Deaths"),
    y = "Cumulative Fatalities"
  )

ggplotly(deathPlot)

As with the spread of the virus, fatalities also grow at an exponential rate. Any slowing in the growth of fatalities is an accurate marker for when the spread of the virus has definitively slowed, despite being a significantly lagging marker. Cumulative fatalities are only shown for countries which have seen 10 or more fatalities for a period beyond 9 days, and for countries with a population > 5,000,000. The blue dashed lines indicate the current position of Australia on this plot.

Daily Increase in Fatalities Vs Confirmed Fatalities

minCases <- 10
minPop <- 10e6
minDays <- 7
dailyVsTotalDeaths <- deaths %>%
  group_by(Country, date) %>%
  summarise_at(vars(deaths), sum) %>%
  dplyr::filter(deaths > 0) %>%
  split(f = .$Country) %>%
  lapply(function(x, n = 7){
    x %>%
      mutate(
        d = c(0, diff(deaths))
      ) %>%
      mutate_at(
        vars(deaths, d), sma, n = n
      )
  }
  ) %>%
  bind_rows() %>%
  inner_join(wikiPops) %>%
  dplyr::filter(deaths > minCases) %>%
  mutate_at(vars(deaths, d), round, 2) %>%
  mutate(n = n()) %>%
  dplyr::filter(
    n > minDays,
    Population > minPop,
    d > 0
  ) %>%
  rename_all(str_to_title) %>%
  rename(
    `Average Daily Fatalities` = D,
    `Average Total Fatalities` = Deaths
  ) %>%
  arrange(Date) %>%
  ungroup() %>%
  mutate(
    Population = comma(round(Population, -3)),
    Country = fct_inorder(Country)
    ) %>%
  ggplot(aes(`Average Total Fatalities`, `Average Daily Fatalities`, colour = Country, label = Date, key = Population)) +
  geom_hline(
    aes(yintercept = `Average Daily Fatalities`),
    data = . %>%
      dplyr::filter(Country == "Australia", Date == max(Date)),
    linetype = 3,
    colour = "blue",
    size = 1/3
  ) +
  geom_vline(
    aes(xintercept = `Average Total Fatalities`),
    data = . %>%
      dplyr::filter(Country == "Australia", Date == max(Date)),
    linetype = 3,
    colour = "blue",
    size = 1/3
  ) +
  geom_point() +
  geom_line() +
  scale_x_log10(labels = comma) +
  scale_y_log10(labels = comma) +
  labs(
    x = "Total Confirmed Fatalities",
    y = 'Daily Fatalities'
  )

capDailyVsTotal <- glue(
  "*Daily Fatalities plotted against Total Fatalities.
  These two values are directly proportional until interventions are successful, at which point the proportional relationship changes, as evidenced by a sudden downwards turn.
  Values shown are 7-day moving averages in order to minimise the impact of day-to-day variation.
  Countries are only shown from the point at which the moving average exceeds {minCases} cases, and have exceeded this value for > {minDays} days.
  Data is additionally restricted to countries with a population > {comma(minPop)}, for greater clarity.
  Australia's position is indicated by the intersection of the dashed blue lines.
  Importantly, __due to the time taken from the initial infection to the day of death, this is a lag indicator of the control of infection__.*
  "
)

ggplotly(dailyVsTotalDeaths)

Daily Fatalities plotted against Total Fatalities. These two values are directly proportional until interventions are successful, at which point the proportional relationship changes, as evidenced by a sudden downwards turn. Values shown are 7-day moving averages in order to minimise the impact of day-to-day variation. Countries are only shown from the point at which the moving average exceeds 10 cases, and have exceeded this value for > 7 days. Data is additionally restricted to countries with a population > 10,000,000, for greater clarity. Australia’s position is indicated by the intersection of the dashed blue lines. Importantly, due to the time taken from the initial infection to the day of death, this is a lag indicator of the control of infection.

Australian Specific Data

ausPops <- tribble(
  ~State, ~Population,
  "New South Wales",    8117976,
  "Victoria", 6629870,
   "Queensland", 5115451,
  "South Australia", 1756494,
  "Western Australia", 2630557,
  "Tasmania", 535500,
  "Northern Territory", 245562,
  "Australian Capital Territory", 428060
)

Australian State populations were taken from the ABS Website and were accurate in Sept 2019. The difference with previous estimates used above was within 0.04%, and as such no adjustments were made.

A series of complimentary charts regarding the spread of COVID-19 are available from the ABC website.

Australia’s spread of the virus appears in the previous plots as marginally slower than many other countries, grouping together with countries such as Greece, Israel and the UK.

• This may be a reflection of Australian geography or our tendency to have large personal space requirements
• This may also be a reflection of the pre-emptive strategies already being taken by the population and government
• Despite the above, Australia is still in an exponential growth phase, however there are early signs of the growth levelling out. The contribution of the distortions caused by the Ruby Princess to this apparent levelling our is still unclear
• Using an estimated population size of 25,203,198, the total percentage of the Australian population confirmed as infected currently sits at 0.022%, or one person in every 4,447.

confirmed %>% 
  dplyr::filter(
    Country == "Australia", 
    date >= max(date) - 1
  ) %>%
  left_join(
    latestAU %>%
      dplyr::select(`Province/State`, date, deaths, recovered)
  ) %>%
  bind_rows(
    group_by(., date) %>%
      summarise_at(vars(confirmed, deaths, recovered), sum, na.rm = TRUE) %>%
      mutate(`Province/State` = "Total")
  ) %>%
  mutate(
    FRate = deaths / confirmed,
    RRate = recovered / confirmed
  ) %>% 
  group_by(`Province/State`) %>% 
  mutate(
    d = diff(confirmed),
    p = diff(confirmed) / min(confirmed)
  ) %>%
  pivot_wider(
    id_cols = "Province/State", 
    names_from = date, 
    values_from = c(confirmed, d, p, deaths, recovered, FRate, RRate)
  ) %>% 
  set_names(
    str_remove(names(.), "confirmed_")
  ) %>%
  dplyr::select(
    1:4, seq(7, ncol(.), by = 2)
  ) %>%
  set_names(
    str_replace(names(.), "^d_.+", "Daily Change")
  ) %>%
  set_names(
    str_replace(names(.), "p_.+", "Percent")
  ) %>%
  set_names(
    str_replace(names(.), "deaths.+", "Fatalities")
  ) %>%
  set_names(
    str_replace(names(.), "recov.+", "Recovered")
  ) %>%
  set_names(
    str_replace(names(.), "FRate.+", "Fatality Rate")
  ) %>%
  set_names(
    str_replace(names(.), "RRate.+", "Recovery Rate")
  ) %>%
  arrange(desc(Percent)) %>%
  ungroup() %>%
  mutate(isTotal = grepl("Total", `Province/State`)) %>%
  arrange(isTotal, desc(Percent)) %>%
  dplyr::select(-isTotal) %>%
  mutate_at(
    vars(Percent, `Fatality Rate`, `Recovery Rate`),
    percent,
    accuracy = 0.1
  ) %>%
  rename(`% Change` = Percent) %>%
  dplyr::select(1:5, starts_with("Fatal"), starts_with("Recov")) %>%
  pander(
    justify = "lrrrrrrrr",
    caption = paste(
      "*Confirmed cases, fatalities and recoveries reported by each state at time of preparation.",
      "Some states are not reporting recovered patient numbers, and these values are presented as empty cells.",
      "Any states with unchanged, or decreasing confirmed cases may indicate delays with the automated data sources, such as health.gov.au or JHU, or that these states have not yet reported for the day.*"
    ),
    emphasize.strong.rows = nrow(.)
  )

Confirmed cases, fatalities and recoveries reported by each state at time of preparation. Some states are not reporting recovered patient numbers, and these values are presented as empty cells. Any states with unchanged, or decreasing confirmed cases may indicate delays with the automated data sources, such as health.gov.au or JHU, or that these states have not yet reported for the day.

Province/State	2020-04-04	2020-04-05	Daily Change	% Change	Fatalities	Fatality Rate	Recovered	Recovery Rate
New South Wales	2,493	2,580	87	3.5%	16	0.6%
Tasmania	80	82	2	2.5%	2	2.4%	15	18.3%
Victoria	1,115	1,135	20	1.8%	9	0.8%	573	50.5%
Queensland	900	909	9	1.0%	3	0.3%
Australian Capital Territory	93	93	0	0.0%	2	2.2%	18	19.4%
Northern Territory	26	26	0	0.0%	0	0.0%	1	3.8%
South Australia	407	407	0	0.0%	0	0.0%	46	11.3%
Western Australia	436	436	0	0.0%	3	0.7%	92	21.1%
Total	5,550	5,668	118	2.1%	35	0.6%	745	13.1%

Plot of Current Australian Values

ausStatsPlot <- ggplotly(
  confirmed %>%
    dplyr::filter(Country == "Australia", confirmed > 0) %>%
    left_join(deaths) %>%
    left_join(recovered) %>%
    group_by(Country, date) %>%
    summarise_at(vars(confirmed, deaths, recovered), sum) %>%
    ungroup() %>%
    mutate(active = confirmed - deaths - recovered) %>%
    pivot_longer(
      cols = c(active, confirmed, deaths, recovered),
      names_to = "Type",
      values_to = "Total"
    ) %>%
    arrange(Type, date) %>%
    split(f = .$Type) %>%
    lapply(mutate, keep = c(TRUE, diff(Total) != 0)) %>%
    bind_rows() %>%
    dplyr::filter(keep, date > ymd("2020-02-28")) %>%
    dplyr::select(-keep) %>%
    mutate(
      Type = str_to_title(Type),
      Type = str_replace_all(Type, "Deaths", "Fatalities"),
    ) %>%
    dplyr::filter(Total > 0) %>%
    rename_all(str_to_title) %>%
    ggplot(aes(Date, Total, colour = Type)) +
    geom_point() +
    geom_smooth(method = "loess", se = FALSE, size = 1/2) +
    scale_y_log10() +
    scale_colour_manual(
      values = c(
        Active = rgb(0, 0, 0),
        Confirmed = rgb(0, 0.3, 0.7),
        Fatalities = rgb(0.8, 0.2, 0.2),
        Recovered = rgb(0.2, 0.7, 0.4)
      )
    )
)
ausStatsPlot$x$data[5:8] %<>%
  lapply(function(x){
    x$hoverinfo <- "none"
    x
  })

ausStatsCap <- "*Current confirmed and recovered cases, along with fatalities for Australia only. Active cases are shown as confirmed cases excluding fatalities and those classed as recovered. Loess curves through all points are shown as continuous lines. Data is only shown from 29^th^ Feb 2020 as data is noticeably sparse prior to this date.*"

ausStatsPlot

Current confirmed and recovered cases, along with fatalities for Australia only. Active cases are shown as confirmed cases excluding fatalities and those classed as recovered. Loess curves through all points are shown as continuous lines. Data is only shown from 29^th Feb 2020 as data is noticeably sparse prior to this date.

State-wise Comparison of Confirmed Infection Rates

minRate <- 5
p5 <- confirmed %>%
  dplyr::filter(
    Country == "Australia"
  ) %>%
  dplyr::rename(State = `Province/State`) %>%
  left_join(ausPops) %>%
  mutate(
    Rate = round(1e6*confirmed / Population, 2),
    Date = format.Date(date, "%d-%B")
  ) %>%
  dplyr::filter(
    !is.na(Population),
    Rate > minRate
  ) %>%
  arrange(date) %>%
  mutate(
    State = fct_inorder(State)
  ) %>%
  dplyr::rename(
    Confirmed = confirmed
  ) %>%
  distinct(State, Confirmed, .keep_all = TRUE) %>%
  ggplot(
    aes(
      x = date, y = Rate, colour = State, 
      label = Date, key = Confirmed
      )
  ) +
  geom_point() +
  geom_smooth(
    method = "lm",
    se = FALSE,
    show.legend = FALSE,
    size = 1/2
  ) +
  geom_line(linetype = 3) +
  scale_y_log10() +
  labs(
    x = "Date",
    y = "Confirmed Infection Rate (cases/million)"
  )
# Hide the tooltip from the regression lines
n <- length(levels(p5$data$State))
p5 <- ggplotly(p5, tooltip = c("Date", "Rate", "State", "Confirmed"))
regIndex <- seq(n + 2, length.out = n, by = 1)
p5$x$data[regIndex] <- lapply(
  p5$x$data[regIndex],
  function(x){
  x$hoverinfo <- "none"
  x
})

p5

Infection rates for each state with data beginning for each state once 5 confirmed cases /million was exceeded. Linear regression lines are shown for each state as solid lines, with NSW perhaps showing a slightly increased rate of infection within the population. Once again, the y-axis is on a logarithmic scale indicating exponential growth is occurring. States are shown in order of the initial date they passed 5 cases/million. Most states show a strong upward curve during the week following 19^th March, clearly showing the impact of the Ruby Princess. The hypothetical growth in the infection rate without that event occurring is not calculated here, however it is clear that the spread of COVID-19 has not continued at the same pace after that event initially occured.

Statistical Analysis of State-wise Infection Rates

lm <- confirmed %>%
  dplyr::filter(
    Country == "Australia"
  ) %>%
  dplyr::rename(State = `Province/State`) %>%
  left_join(ausPops) %>%
  mutate(
    Rate = round(1e6*confirmed / Population, 2),
    Date = format.Date(date, "%d-%B")
  ) %>%
  dplyr::filter(
    !is.na(Population),
    Rate > minRate,
    date <= dt
  ) %>%
  arrange(desc(confirmed)) %>%
  mutate(
    State = fct_inorder(State)
  ) %>%
  dplyr::rename(
    Confirmed = confirmed
  ) %>%
  arrange(State, Date) %>%
  distinct(State, Confirmed, .keep_all = TRUE) %>%
  with(
    lm(log10(Rate) ~ (State + date)^2)
  ) 

In order to test whether the infection rates are different between states, a linear regression model was fit. \(\log_{10}\)(Cumulative Confirmed Infection Rate) was assigned as the response variable with predictor variables being the State and Date. Each State was assigned its own intercept and slope by use of an interaction term (i.e. State:date). Given the potentially larger slope in NSW, this state was set as the baseline, with each other slope (i.e. interaction term) being presented as the difference in slope between each state and NSW. In this way comparisons against NSW were performed, but no comparisons between other states were performed.

Differences in the State-level intercepts are not particularly relevant, apart from capturing the initial time at which cumulative confirmed infection rate exceeded 5 cases / million. Differences between State-level slopes and NSW however, are of great interest, and as such only the slopes are shown. For NSW this captures the actual slope of the daily change in infection rate, whilst for all other States, this represents the difference between the daily change in infection rate for that State in comparison to NSW. Only differences in slope with an Bonferroni-adjusted p-value < 0.05 are of particular interest. For those States which appear to be of interest, a negative value indicates an infection rate increasing more slowly than NSW, whilst a positive value indicates the opposite. If no difference is noted between any state and NSW, it can be consideed that these states are all experiencing growth at a similar rate. Currently, the infection rate in most states is doubling every 4.2 days.

st <- shapiro.test(resid(lm))$p.value
bp <- olsrr::ols_test_breusch_pagan(lm)$p
lm %>%
  tidy() %>%
  mutate(
    term = str_remove_all(term, "State")
  ) %>%
  rename(
    Term = term,
    Estimate = estimate,
    SE = std.error,
    `T` = statistic,
    p = p.value
  ) %>%
  dplyr::filter(
    str_detect(Term, "date")
  ) %>%
  mutate(
    Term = str_remove(Term, ":date"),
    Term = str_replace_all(Term, "date", "NSW Daily Increase"),
    adjP = p.adjust(p, method = "bonf"),
    Estimate = sprintf("%.4f", Estimate),
    SE = sprintf("%.3f", SE),
    `T` = sprintf("%.2f", `T`),
    p = case_when(
      p < 1e-4 ~ sprintf("%.2e", p),
      p >= 1e-4 ~ sprintf("%.4f", p)
    ),
    adjP = case_when(
      adjP < 1e-4 ~ sprintf("%.2e", adjP),
      adjP >= 1e-4 ~ sprintf("%.4f", adjP)
    )
  ) %>%
  rename(
    `p~bonf~` = adjP
  ) %>%
  pander(
    justify = "lrrrrr",
    emphasize.strong.rows = which(
      as.numeric(.$`p~bonf~`) < 0.1 & .$Term != "date"
    ),
    caption = paste(
      "*Results of linear regression analysis",
      "comparing the slopes of lines which track __daily",
      "change in cumulative confirmed infection rates__,",
      "scaled for population size in each state.",
      "Intercept terms are not shown.",
      "The daily slope of the NSW regression line is given as the first value, with any differences between this line and other states shown as the subsequent rows.",
      "All slopes and differences in slopes are given on the log~10~ scale.",
      "Any highlighted state indicates an infection rate which is increasing at a different daily rate to NSW.",
      "Positive values in the Estimate column indicate a greater rate the NSW, whilst a negative value in the Estimate column indicates a slower growth rate than NSW.",
      "The fitted model was also checked using the Shapiro-Wilk Test for normality",
      paste0(
        "(p = ",
        sprintf( "%.4f", st),
        "). ",
        ifelse(st < 0.05, "The current results indicate potential violations of the assumption of normality and results in this table __should be taken with caution__.", "No issues with the assumption of normality were noted.")
      ),
      paste0(
        "Heteroscdasticity was tested using the Breusch-Pagan test (p = ",
        sprintf( "%.4f", bp),
        "). ",
        ifelse(bp < 0.05, "Given that unequal variances were detected, results and predictions from this model __should be taken with caution__* The source of this unequal variance is the spike in cases immediately following the Ruby Princess disembarkation.", "No issues with the assumption of constant variance were noted.*")
      )
    )
  )

Results of linear regression analysis comparing the slopes of lines which track daily change in cumulative confirmed infection rates, scaled for population size in each state. Intercept terms are not shown. The daily slope of the NSW regression line is given as the first value, with any differences between this line and other states shown as the subsequent rows. All slopes and differences in slopes are given on the log₁₀ scale. Any highlighted state indicates an infection rate which is increasing at a different daily rate to NSW. Positive values in the Estimate column indicate a greater rate the NSW, whilst a negative value in the Estimate column indicates a slower growth rate than NSW. The fitted model was also checked using the Shapiro-Wilk Test for normality (p = 0.2954). No issues with the assumption of normality were noted. Heteroscdasticity was tested using the Breusch-Pagan test (p = 0.0636). No issues with the assumption of constant variance were noted.

Term	Estimate	SE	T	p	pbonf
NSW Daily Increase	0.0721	0.003	21.58	1.95e-45	1.56e-44
Victoria	-0.0026	0.005	-0.50	0.6207	1.0000
Queensland	-0.0065	0.005	-1.24	0.2155	1.0000
Western Australia	0.0005	0.006	0.10	0.9242	1.0000
South Australia	0.0012	0.005	0.23	0.8180	1.0000
Australian Capital Territory	0.0133	0.007	1.98	0.0493	0.3941
Tasmania	-0.0108	0.005	-1.99	0.0489	0.3913
Northern Territory	-0.0139	0.009	-1.60	0.1127	0.9015

new <- confirmed %>% 
  dplyr::filter(Country == "Australia", date == dt) %>%
  mutate(date = dt + 1) %>%
  dplyr::select(State = `Province/State`, date, confirmed) %>%
  droplevels()
new %>%
  cbind(
    predict.lm(
      lm,
      newdata = new, 
      interval = "predict",
      level = 0.95
    )
  ) %>%
  left_join(ausPops) %>%
  pivot_longer(
    cols = c(fit, lwr, upr),
    names_to = "type",
    values_to = "value"
  ) %>%
  mutate(
    value = round(Population * 10^value / 1e6, 0)
  ) %>%
  pivot_wider(
    names_from = type,
    values_from = value
  ) %>%
  arrange(desc(Population)) %>%
  mutate(
    lwr = case_when(
      lwr < confirmed ~ confirmed,
      lwr >= confirmed ~ lwr
    )
  ) %>%
  mutate(PI = glue("[{lwr}, {upr}]")) %>%
  dplyr::select(State, Confirmed = confirmed, PI) %>%
  set_names(
    case_when(
      names(.) == "Confirmed" ~ as.character(dt),
      names(.) == "PI" ~ paste0("95% Pred Interveral (", as.character(dt + 1), ")"),
      names(.) != "Confirmed" ~ names(.)
    )
  ) %>%
  pander(
    justify = "lrr",
    caption = glue(
      "*95% prediction intervals for cumulative confirmed cases in each state on {format(dt, '%B %d, %Y')}.
      Values were obtained using the above linear model, and as such, any caveats about the appropriateness of the model should be taken into account.*"
    )
  )

95% prediction intervals for cumulative confirmed cases in each state on April 04, 2020. Values were obtained using the above linear model, and as such, any caveats about the appropriateness of the model should be taken into account.

State	2020-04-04	95% Pred Interveral (2020-04-05)
New South Wales	2,493	[2589, 9118]
Victoria	1,115	[1115, 3857]
Queensland	900	[900, 3163]
Western Australia	436	[463, 1682]
South Australia	407	[444, 1588]
Tasmania	80	[80, 256]
Australian Capital Territory	93	[113, 412]
Northern Territory	26	[26, 64]

Assessment of Testing Within Each State

tested <- list.files("tested", pattern = "tsv", full.names = TRUE) %>%
  sort() %>%
  .[length(.)] %>%
  read_tsv()

tested <- c(
  "New South Wales" = nswTests,
  "Queensland" = qldTests,
  "Victoria" = vicTests,
  "Western Australia" = waTests,
  "South Australia" = saTests,
  "Australian Capital Territory" = actTests,
  "Tasmania" = tasTests,
  "Northern Territory" = ntTests
) %>%
  enframe(
    name = "State", value = "Tested"
  ) %>%
  mutate(Tested = as.numeric(Tested)) %>%
  dplyr::filter(!is.na(Tested)) %>%
  bind_rows(tested) %>%
  group_by(State) %>%
  summarise(Tested = max(Tested, na.rm = TRUE))
tested %>%
  write_tsv(
    glue(
      "tested/tested_{format(Sys.Date(), '%Y%m%d')}.tsv"
    )
  )

Testing numbers were initially sourced from manual inspection of individual press releases for NSW, QLD, VIC, WA, SA, TAS, ACT and the NT. Updates on testing numbers beyond the initial values were performed automatically using the above code which probes each state’s official releases for the latest values, and updates where found. The number of tested individuals in each state was then assessed as a function of State population size. All results are valid at the time of report generation.

tested %>%
  full_join(
    confirmed %>%
      dplyr::filter(
        date == max(date)
      ) %>%
      rename(
        State = `Province/State`
      )
  ) %>%
  mutate(
    Tested = case_when(
      is.na(Tested) ~ confirmed,
      Tested < confirmed ~ confirmed,
      !is.na(Tested) ~ Tested
    )
  ) %>%
  left_join(ausPops) %>%
  bind_rows(
    tibble(
      State = "Total",
      Population = sum(.$Population, na.rm = TRUE),
      confirmed = sum(.$confirmed, na.rm = TRUE),
      Tested = sum(.$Tested, na.rm = TRUE)
    )
  ) %>%
  mutate(
    `Percent Tested` = Tested / Population,
    Positive = confirmed / Tested,
    Negative = 1 - Positive,
    isTotal = grepl("Total", State)
  ) %>%
  dplyr::select(
    State, Population,
    Confirmed = confirmed,
    Tested, 
    contains("Percent"), 
    ends_with("ive"),
    isTotal
  ) %>%
  arrange(isTotal, desc(`Percent Tested`)) %>%
  dplyr::select(-isTotal) %>%
  dplyr::rename(
    `% Tested` = `Percent Tested`,
    `% Positive Tests` = Positive,
    `% Negative Tests` = Negative
  ) %>%
  mutate_at(
    vars(starts_with("%")), percent, accuracy = 0.01
  ) %>%
  pander(
    justify = "lrrrrrr",
    missing = "",
    caption = glue(
      "*COVID-19 testing scaled by state population size.
      Confirmed cases are assumed to be the tests returning a positive result.
      The current numbers available for some states are a lower limit, and as such, the proportion of the population tested is likely to be higher, as is the proportion of tests returning a negative result.*"
    ),
    emphasize.strong.rows = nrow(.)
  )

COVID-19 testing scaled by state population size. Confirmed cases are assumed to be the tests returning a positive result. The current numbers available for some states are a lower limit, and as such, the proportion of the population tested is likely to be higher, as is the proportion of tests returning a negative result.

State	Population	Confirmed	Tested	% Tested	% Positive Tests	% Negative Tests
South Australia	1,756,494	407	32,000	1.82%	1.27%	98.73%
New South Wales	8,117,976	2,580	117,633	1.45%	2.19%	97.81%
Australian Capital Territory	428,060	93	5,328	1.24%	1.75%	98.25%
Northern Territory	245,562	26	2,642	1.08%	0.98%	99.02%
Queensland	5,115,451	909	51,785	1.01%	1.76%	98.24%
Victoria	6,629,870	1,135	54,000	0.81%	2.10%	97.90%
Western Australia	2,630,557	436	17,251	0.66%	2.53%	97.47%
Tasmania	535,500	82	1,779	0.33%	4.61%	95.39%
Total	25,459,470	5,668	282,418	1.11%	2.01%	97.99%

R Session Information

R version 3.6.3 (2020-02-29)

Platform: x86_64-pc-linux-gnu (64-bit)

locale: LC_CTYPE=en_AU.UTF-8, LC_NUMERIC=C, LC_TIME=en_AU.UTF-8, LC_COLLATE=en_AU.UTF-8, LC_MONETARY=en_AU.UTF-8, LC_MESSAGES=en_AU.UTF-8, LC_PAPER=en_AU.UTF-8, LC_NAME=C, LC_ADDRESS=C, LC_TELEPHONE=C, LC_MEASUREMENT=en_AU.UTF-8 and LC_IDENTIFICATION=C

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: plotly(v.4.9.2), DT(v.0.12), pander(v.0.6.3), RCurl(v.1.98-1.1), rvest(v.0.3.5), xml2(v.1.2.2), jsonlite(v.1.6.1), glue(v.1.3.1), broom(v.0.5.4), ggrepel(v.0.8.1), matrixStats(v.0.55.0), scales(v.1.1.0), lubridate(v.1.7.4), magrittr(v.1.5), forcats(v.0.4.0), stringr(v.1.4.0), dplyr(v.0.8.4), purrr(v.0.3.3), readr(v.1.3.1), tidyr(v.1.0.2), tibble(v.2.1.3), ggplot2(v.3.2.1) and tidyverse(v.1.3.0)

loaded via a namespace (and not attached): nlme(v.3.1-144), bitops(v.1.0-6), fs(v.1.3.1), httr(v.1.4.1), rprojroot(v.1.3-2), tools(v.3.6.3), backports(v.1.1.5), R6(v.2.4.1), nortest(v.1.0-4), DBI(v.1.1.0), lazyeval(v.0.2.2), colorspace(v.1.4-1), withr(v.2.1.2), gridExtra(v.2.3), tidyselect(v.1.0.0), curl(v.4.3), compiler(v.3.6.3), cli(v.2.0.1), Cairo(v.1.5-10), labeling(v.0.3), goftest(v.1.2-2), digest(v.0.6.25), foreign(v.0.8-75), rmarkdown(v.2.1), rio(v.0.5.16), pkgconfig(v.2.0.3), htmltools(v.0.4.0), dbplyr(v.1.4.2), fastmap(v.1.0.1), highr(v.0.8), htmlwidgets(v.1.5.1), rlang(v.0.4.4), readxl(v.1.3.1), rstudioapi(v.0.11), shiny(v.1.4.0), generics(v.0.0.2), farver(v.2.0.3), zoo(v.1.8-7), crosstalk(v.1.0.0), zip(v.2.0.4), car(v.3.0-6), Rcpp(v.1.0.3), munsell(v.0.5.0), fansi(v.0.4.1), abind(v.1.4-5), lifecycle(v.0.1.0), stringi(v.1.4.6), yaml(v.2.2.1), carData(v.3.0-3), MASS(v.7.3-51.5), grid(v.3.6.3), promises(v.1.1.0), crayon(v.1.3.4), lattice(v.0.20-40), haven(v.2.2.0), hms(v.0.5.3), knitr(v.1.28), pillar(v.1.4.3), reprex(v.0.3.0), evaluate(v.0.14), data.table(v.1.12.8), modelr(v.0.1.6), olsrr(v.0.5.3), vctrs(v.0.2.3), httpuv(v.1.5.2), selectr(v.0.4-2), cellranger(v.1.1.0), gtable(v.0.3.0), assertthat(v.0.2.1), openxlsx(v.4.1.4), xfun(v.0.12), mime(v.0.9), xtable(v.1.8-4), later(v.1.0.0), viridisLite(v.0.3.0), ellipsis(v.0.3.0) and here(v.0.1)

Report prepared by the University of Adelaide, Bioinformatics Hub. https://www.adelaide.edu.au/bioinformatics-hub/